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ABSTRACT 

We compare the statistics of driven, supersonic turbulence at high Mach number using FLASH 
a widely used Eulerian grid-based code and PHANTOM, a Lagrangian smoothed particle hy- 
drodynamics (SPH) code at resolutions of up to 512 3 in both grid cells and SPH particles. We 
find excellent agreement between codes on the basic statistical properties: a slope of fc -195 in 
the velocity power spectrum for hydrodynamic, Mach 10 turbulence, evidence in both codes 
for a K olmogorov-like slope of fc~ 5 / 3 in the variable p 1 / 3 v as suggested by Krits uk et alj 
(2007|) and a log-normal PDF with a width that scales with Mach number and proportion- 
ality constant b = 0.33 — 0.5 in the density variance-Mach number relation. The measured 
structure function slopes are not converged in either code at 512 3 elements. 

We find that, for measuring volumetric statistics such as the power spectrum slope and 
structure function scaling, SPH and grid codes give roughly comparable results when the 
number of SPH particles is approximately equal to the number of grid cells. In particular, to 
accurately measure the power spectrum slope in the inertial range, in the absence of sub-grid 
turbulence models, requires at least 512 3 computational elements in either code. On the other 
hand the SPH code was found to be better at resolving dense structures, giving maximum 
densities at a resolution of 128 3 particles that were similar to the maximum densities resolved 
in the grid code at 512 3 cells, reflected also in the high density tail of the PDF. We find SPH 
to be more dissipative at comparable numbers of computational elements in statistics of the 
velocity field, but correspondingly less dissipative than the grid code in the statistics of density 
weighted quantities such as p 1</3 v. 

For SPH simulations of high Mach number turbulence we find it important to use suffi- 
cient non-linear /3-viscosity in order to prevent particle interpenetration in shocks (we require 
Pvisc = 4 instead of the widely used default value, f3 V i SC = 2). 

Key words: hydrodynamics - Interstellar Medium (ISM) - methods: numerical - shock 
waves - stars: formation - turbulence 



1 INTRODUCTION 

Dense interstellar molecular clouds are ubiquitously observed to 
have non-thermal line widths im plying supersonic internal mo- 
tions dZuckerman & Evans! 1 19741) . Furthermore the amplitude of 
such motions increases with spatial s cale in a manner reminiscent 
of turbulent flows in the lab oratory ((Larson 198ll: ISolomon et all 
1 19871 : iHever & Bruntl |2004|) . Understanding the nature and ori- 
gin of such 'supersonic turbulence' is theref ore key — perhaps 
the key — to understanding star formation (jEl megreen & Scald 
120041 ; [Mac Low & Klessenl |2004| ; iMcKee & Ostrikerl l2007h . Tur- 
bulence provides a natural explan ation for the clustered and hier - 
archical nature of star formation dElmegreen & Falgarone 19961): 
the measured fractal dimen sion of interstellar gas dKritsuk et alj 
120071 : iFederrath et alj|2009l and references therein); the few per- 



cent e fficiency with which gas is converted into stars iPadoanl 



19951: Vazquez-Semadeni. Ballesteros-Paredes & Klessen 20031 : 



Krumholz & McKed 120051: lElmegreenl 120081) : most likely deter- 



mines in large part the mass dis tribution of star forming cores 
(the core mass d istribution, CMD) dBallesteros-Paredes et al . 2006; 
iDib et alj [2008) and possibly the mass distrib ution of stars them- 
selves, i.e., the Initial Ma s s Function (IMF) l Padoan et~al] 1 1 9971 : 
IPadoan & Nordlundl |2002| ; iHennebelle & Chabrien 1 |2008h . How- 
ever, given that a full theory of turbulence is elusive even in the 
incompressible reg ime apart from the phenomenology provided by 
iKolmogorovl dl94ll) . one inevitably turns to numerical simulations 
to glean insight. 

Given the importance of numerical simulations in understand- 
ing the basic statistics of supersonic turbulence, and the possible 
implications for star formation theory, it is crucial that results in- 
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ferred from such simulations are robust with respect to different 
numerical methods and codes. This has motivated at least two ma- 
jor code comparison projects in the last year or so in which both of 
the present authors ha ve been involved. The 'Potsdam comparison' 
dKitsionas et ail [2009) compared simulations of decaying, hydro- 
dynamic turbulence using 7 different codes (3 SPH codes and 4 
grid-based codes) at a fixed resolution (215 3 particles for the SPH 
codes, and 256 3 grid cells for the grid codes). The results showed 
generally good agreement on statistics suc h as the density PDFs 
and power spectra, similar to earlier studies ( Ma c Low et al.|[l998|: 
lKlessenetal1l2000h . In a similar spirit the KITP07 comparisorU 
compared a large number of grid based and SPH calculations of 
decaying turbulence, for both hydrodynamics and MHD, at a range 
of resolutions, though the results are yet to be published. 

Both of these comparisons are problematic in several respects. 
The first is that it is difficult to make a statistical comparison us- 
ing decaying turbulence, since the time evolution is limited and 
therefore only a few instantaneous snapshots can be compared. 
Instantaneous snapshots however are subject to intermittent fluc- 
tuations that make a head-to-head comparison based on single 
time slices difficult dKritsuk et alj2007l : lFederrafh et al.l2010h . Sec- 
ondly, both comparisons start from evolved initial conditions, pro- 
duced from either a previously driven SPH simulation (Potsdam) or 
a grid-based calculation (KITP) that has to be interpolated and/or 
downsampled to/from the grid/particles appropriate to the differ- 
ent codes, with an ensuing loss of accuracy and consistency before 
the comparison has even begun (this problem is much worse for the 
MHD case where differences in divergence-free representations for 
the magnetic field between codes is a further issue). 

In this paper, we consider only two codes, an SPH code, 
PHANTOM, and a grid-based code, FLASH, which we take to be 
broadly representative of the fundamentally different classes of 
code used for star formation studies (the codes are described in 
fQ. The turbulence in both codes is driven from stationary, uniform 
initial conditions with exactly the same energy input and driving 
pattern over multiple turbulent crossing times. We also consider a 
range of resolutions (128 3 , 256 3 and 512 3 in both grid elements 
and the number of SPH particles) in order to estimate resolution 
requirements and establish where convergence has occurred in one 
code or the other, or neither. In the present work we limit ourselves 
to a study of hydrodynamic turbulence, that is, without magnetic 
fields. This is primarily because the algorithms for Magnetohydro- 
dyn amics in SPH curre ntly being used for star formation studies 
(e.g. |Price&Batel2008h rely on the Euler potentials formulation of 
the magnetic field, that cannot be used for turbulence studies over 
multip le crossin g times due to the restricted field representations 
(see |Pricell2010l for recent progress). 

The goals of this paper are to: i) establish whether or not agree- 
ment can be found between SPH and grid codes on the basic statis- 
tics of supersonic turbulence; ii) define resolution criteria for var- 
ious statistical measures of supersonic turbulence such as power 
spectra, PDFs and structure functions; and iii) establish the relative 
strengths and weaknesses of each method for turbulence studies. 
We discuss the numerical methods in S|2] with the Fourier space 
driving discussed in $2.2\ Results from both codes are presented in 
!j3]and our findings discussed in [j4] 



2 NUMERICAL METHOD 
2.1 Equations 

We solve the equations of non-self-gravitating hydrodynamics 
given by 



| + (v'V)p = -pV-v, 
<9v , „, VP 

— + (V • V V = + f 3 tir 

at p 



(1) 

(2) 



where i s ur is a stirring force, the details of which are discussed 
below ( £|2.21 >- The pressure is related directly to the density via an 
isothermal equation of state 



P = c s p, 



(3) 



where, since the equations are scale-free to all but the Mach num- 
ber, we use c s = 1. We solve using periodic boundary con- 
ditions in the three-dimensional domain x, y, z £ [0, 1]. The initial 
conditions are a uniform density medium p — po — 1 with zero 
initial velocities. 

We discuss our results in terms of the dynamical time, defined 
as td = L /(2M ) , where L is the box size and M is the RMS Mach 
number. However, the absence of a physical scale in the equations 
means that the results can be arbitrarily scaled to the interstellar 
medium by defining length, mass and time scales. For example, 
adopting a length scale of 10 pc and a sound speed of 0.2 km/s, 
gives the physical time from t/td according to 



physical 



2.5 Myr 



10 pc J V0.2 km/s 



-i t 



(4) 



One can similarly set a mass scale by defining the initial density 
to be no ~ 3 x 10 2 cm -3 , i.e., po ~ 10~ 21 gem -3 assuming 
fully molecular hydrogen gas, giving a total mass in the box of ~ 
1.5 x 10 4 M©. For these parameters, the maximum density reached 
in our calculations from turbulent fluctuations alone is p ma x ~ 
10~ 17 -10" 16 gem" 3 or n max ~ 10 B cm" 3 (see Fig.|2J for the 
highest resolution SPH calculation. 



2.2 Driving 

F or driving tur b ulence we use the same driving r o utine used 



iilSchmidt et al 



j 200i) : lFederrathetal.l J2008L l201Cl l2009h and 
2009h . The driving routine updates a vector of 



ISchmidt etail 

real values according to an algorithm that gener ates an Ornstein- 
Uhlen beck, or "coloured noise" sequence (e.g. Eswaran & Pope 
1988). The sequence x n is a Markov process that takes the previ- 
ous value, weights by an exponential damping factor with a given 
auto-correlation time t s , and drives by adding a Gaussian random 
variable, weighted by a second damping factor, also with correla- 
tion time t s . For a timestep dt, this sequence can be written as: 



-l = fx„ + oV (1 - f 2 ) z n 



(5) 



1 http://kitpstarformation07.wikispaces. com/Star +Formation+Test+Problems 



where / = exp(— dt/t s ), and z n is a Gaussian random variable 
drawn from a Gaussian distribution with unit variance, and a is 
the desired var iance of the Ornstein-Uhlenbeck sequence (see e.g. 
lBartoscrJ|200lh . The resulting sequence satisfies the properties of 
zero mean, and stationary RMS equal to a. Its power spectrum in 
the time domain can vary from white noise (P(f) = const) to 
"brown" noise (P(f) oc l// 2 ). 

The physical forcefield is constructed in Fourier space using 
the Ornstein-Uhlenbeck process. This allows a simple decomposi- 
tion of the field into a solenoidal (divergence-free) part and a com- 
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pressible (curl-free) part using a Helmholtz decomposition. In this 
study, we only keep the solenoidal part. Inverse Fourier transforma- 
tion yields the physical solenoidal force field f s tir used in equation 
l|2}. A more d etailed description of th e forcing module applied here 
is provided in iFederrath et allbOlOl) . 

The time-dependent Fourier modes for constructing the forc- 
ing patterns f s u r were calculated and written to a file before the 
actual numerical experiments. Both the SPH and the grid code read 
exactly the same forcing sequence from this file. Thus, it was guar- 
anteed that both codes were using exactly the same forcing at all 
times during the comparison experiments. 



2.3 flash (grid) 
2.3.1 Hydrodynamics 

flash dFrvxell et all2000l:lDubev et alj2008h is an adaptive-mesh 
refinement code dBerger & Colella 1989) that uses t he piecewise 
parabolic method (PPM. IColella & Woodwardl 19841) to solve the 
equations of hydrodynamics. The PPM provides a shock capturing 
scheme to keep shocks and contact discontinuities sharp (typically 
spreading over 2-3 zones), while maintaining third order accuracy 
in smooth flows through a parabolic reconstruction scheme. In this 
study, FLASH v3 was used, which provides a uniform grid mode. 
Thus, the overhead in storing and iterating the adaptive mesh hier- 
archy was completely removed, which yields a speed-up of factors 
of a few. FLASH is parallelised using the message passing inter- 
face (MPI). For the resolutions studied here (128 3 , 256 3 and 512 3 



grid cells), 1, 8 and 64 MPI processes respectively were used in a 
mode of parallel computation, each calculation taking roughly 12, 
250 and 5000 CPU-hours respectively . FLASH has been extensively 
tested against laboratory experime nts ( Calder et al. 2002) an d other 
codes dPimonte et~al]|2004l : iHeitmann et alj|2005l ; iKitsionas et all 
120091) . 



2.3.2 Tracer particles 

FLASH provides an option for Lagrangian tracer particles, which 
can be evolved alongside the hydrodynamics. Similar to SPH parti- 
cles, tracer particles provide information in the Lagrangian frame, 
but unlike SPH particles, the tracer particles have no feedback on 
the hydrodynamics, i.e. the variables on the grid are independent of 
the tracers. The tracer particles' x, y and z positions can be any real 
number within the computational domain, not bound to the grid. 
However, they are moved with the velocity computed on the grid. 
The velocity is interpolated at the exact position of each tracer par- 
ticle for each timestep using a first order cloud-in-cell interpolation 
scheme. Higher-order interpolation schemes like the triangular- 
shaped-cloud scheme can also be used instead. However, we used 
the first-order scheme here, because various tests suggested no 
strong dependence of our results on the interpolation scheme. The 
tracer particles were moved on the hydrodynamic timestep with the 
grid-interpolated velocity using a first-order scheme. We initialised 
128 3 , 256 3 and 512 3 tracer particles at t = on a uniform grid 
at exactly the same positions as the SPH particles were initialised 
in the PHANTOM calculations (see E|2.4.4b . matching the grid and 
SPH resolutions (128 3 , 256 3 and 512 3 , respectively). Adding the 
tracer particles does not add any significant computational over- 
head to the FLASH calculations, apart from the additional memory 
requirements. 

In order to extract the maximum possible information from 
the tracer particles, we have computed — in post-processing — 



a density field based solely on the tracer particle positions. This 
is achieved by assuming they are particles of fixed mass (divid- 
ing the total mass in the simulation by the number of tracer parti- 
cles) and using the SPH density calculation routine from PHAN- 
TOM where the density and smoothing length are iterated self- 
consistently (based on Eqs.[6]and[7}. Column-integrated and cross- 
section slice plots of the density fi eld were the n produced as for the 
PHANTOM results using SPLASH jPricell2007h . 

2.4 PHANTOM (SPH) 

PHANTOM is a low-memory, highly efficient SPH code written es- 
pecially for studying non-self-gravitating problems. The code is 
made very efficient by using a simple neighbour finding scheme 
based on a fixed grid and linked lists of particles. The calculations 
shown in this paper have used only the shared memory OPENMP 
parallelisation in PHANTOM, using 4, 8 and 32 processors and re- 
quiring 265, 5050 and 120,000 CPU-hours for the 128 3 , 256 3 and 
512 3 calculations, respectively. Thus the 256 3 PHANTOM calcu- 
lation was roughly comparable in computational cost to the 512 3 
FLASH calculation, and similarly for the 128 3 PHANTOM vs. 256 3 
FLASH (though some caution is required here due to the differ- 
ent machines and architectures used to run each code). One may 
also consider that PHANTOM was found to be roughl y an order of 
magni tude faster than 'standard' SPH codes in the lKitsionas et al] 
d2009l) turbulence comparison. 



2.4.1 Hydrodynamics 

For hydrodynamics PHANTOM implements the full vari- 
able smoothing le ngth SPH formulation develo p ed b y 
IPrice & MonagharJ d2004l) and IPrice & Monagharj d2007l) . 
whereby the smoothing length, h, and density, p, are mutually 
dependent via the density sum (for particle a) 



Pa 



y^ y m b Wab{h a ) 



which is an exact solution to 0, and the relation 



h a = V 



Pa 



1/3 



(6) 



(7) 



where m is the particle mass and W a b = iy(|r a — rb\,h a ) is 
the SPH smoothi ng kernel (see e.g. Monaghan 1992: IPricdliool: 
Monaghan 2005 for reviews of SPH). Equations ^ and Q are 
iterated s elf-consistently using a Ne wton-Raphson method as de- 
scribed in lPrice & Monaghan! d2007l) . where in this paper we have 
used rj = 1.2, giving approximately 58 neighbours per particle in a 
smooth distribution. 

The fact that the smoothing length has a functional depen- 
dence on (ultimately) the particle position means that the deriva- 
tives of h can be accounted for in the equations of motion, resulting 
in exact conservation of momentum, angular momentum, energy 
and entropy in the SPH equations. In PHANTOM the equations of 
motion d2} take the form 



+ 



Pg+qa 

Pb + Qb 



V a W ab (h a ) 



bPt 



VaW ab {h b ) 



+ fs; 



(8) 



where P is the pressure, is a dim ensionless quantity related to 
the smoothing length gradients (see IPrice & Monaghan! |2007| for 
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details) and q represents the artificial viscosity term (discussed be- 
low). In the absence of shock dissipation (q = 0) there is zero 
numerical dissipation contained in the above equations and energy 
is conserved to the accuracy of the timestepping scheme — here a 
Kick-Drift-Kick leapfrog integrator equivalent to the velocity Ver- 
let method, implemented with individual particle timesteps. For an 
isothermal equation of state the viscosity term q therefore repre- 
sents the only numerical form of energy loss. 



2.4.2 Artificial viscosity 

Shocks are treated in PHANTOM using a standard artificial viscos- 
ity term, though formulated slightly differently to the usual SPH 
expression in order to obtain a more efficient calculation. Instead 
of the usual expression we write the artificial viscosity as q in (8), 
where 



la = 



where v ab 



■^a a paVsig,ayVab ■ Yab\, 



a — Vf, and we use 



V ai) ■ iab <0 
Va6 ■ f ab > 



,a C-s,a H - /3visc\v ab ' ^*ab\ 



(9) 



(10) 



as the maximum signal velocity for hydrodynamics. The /3 term 
in the signal velocity provides a non-linear term that was origi- 
nally introduced to pr event particle pe netration in high Mach num- 
ber shocks (see e.g. lMonaghan| [l989). Indeed one of our findings 
from this comparison is that sufficient /3-viscosity is an important 
factor for accurate SPH calculations in the supersonic regime. We 
use /3 v isc = 4 in this paper, the motivation for which is discussed 
further in ^3.2. H and demonstrated in Appendix lAl 

The art ificial viscosity d escribed above is essentially the 
same as the iMonaghanl d 19971) formulation with a slightly dif- 
ferent averaging of t he den sity and signal velocity. We use the 
iMorris & Monaghar] d 19971) switch to reduce dissipation away 
from shocks, in which the dissipation parameter a is evolved ac- 
cording to a source and decay equation 



da a 

dt 



+ S a , 



h a /(ac s 



(11) 



where we have used a — 0.1, 5 = max(0, —V ■ v), a m in = 0.05, 
enforced a max ~ 1.0, and given all particles a = a m i n initially. 



2.4.3 Boundary conditions 

Periodic boundary conditions are implemented in PHANTOM by di- 
rectly finding neighbours across the periodic boundary and inter- 
acting with a distance in each direction calculated according to 



(x a 



Xb) = ram 



(x a Xb) ^ 



Xb\ 



L). 



■ x b 



, (12) 



where L is the box size in the corresponding direction. This gives 
a significant memory saving since memory does not have to be as- 
signed to the storage of ghost particles. 



2.4.4 Initial conditions 

The SPH particles were set up initially on a regular cubic lattice, 
using equal mass particles, identical to the setup used for the La- 
grangian tracer particles in the FLASH calculations ( ^2.3.2t . 
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Figure 1. Mass-weighted RMS Mach number as a function of time for the 
six calculations, as indicated in the legend. The time evolution is similar for 
both the SPH and the grid code and for all resolutions up to ~ 44^, where 
all calculations show deviations from each other of order 5% in 8M /M. 
though with no systematic trends with either code or resolution. 
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Figure 2. Maximum density as a function of time in the calculations us- 
ing PHANTOM (SPH, black lines) and FLASH (grid, red lines) at resolutions 
of 128 3 , 256 3 and 512 3 particles/grid cells, as indicated. The evolution 
shows strong time variability, a consequence of the intermittency inherent 
in the log-normal probability distribution function (increasingly higher den- 
sities have correspondingly smaller probabilities). The maximum density is 
therefore also a strong function of resolution in each code. At 128 3 parti- 
cles the SPH code resolves maximum densities similar to those achieved at 
512 3 on the grid. 



3 RESULTS 

Both codes have been run using the same driving pattern for 10 
dynamical times (t = 0.5 in code units) and at a resolution of 
128 3 , 256 3 and 512 3 elements. For the grid code, this resolution 
is fixed spatially throughout the evolution, giving fixed resolution 
in volume but variable resolution in mass, whilst for the SPH code 
the particles move following the fluid motion, giving equal resolu- 
tion in mass but variable resolution in volume. The two methods 
are therefore very nicely complementary for assessing the statis- 
tics of supersonic turbulence for different quantities which may 
be either mass or volume-weighted. It should be noted that whilst 
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grid-based calculations of ISM turbulen ce have been run at much 
higher resolutions — up to 2048 3 , see iKritsuk et alj|2007l with 
resolutio ns of 1024 3 even in early simulations of decaying tur- 
bulence dPorter et alJI 19981) — our use of 134,217,728 SPH par- 
ticles (512 3 ) represents the highest resolution turbulence simula- 
tion performed to date with an SPH code, over an order of mag- 
nitude highe r than the "high resolution" c alculation (10 million 
particles) in Bal lesteros- Paredes et al] fe006l) and two-and-a-half 
orders of magnitude higher than the ~ 200,000 particles used for 
many of the runs in that paper and elsewhere dKlessen et afll200d : 
IVazquez-Semadeni et alj|2003h . 

3.1 Time evolution of global variables 

The time variation of the mass-weighted RMS Mach number and of 
the maximum gas density are shown in Figure s[T] and [2] computed 
at every timestep for both codes. The mass-weighted RMS Mach 
number in SPH is simply the square root of the average value of 
v 2 /c 2 on the particles, whilst in FLASH this has been computed 
using the RMS value of pv 2 /(pod)), to correspond to the SPH av- 
erage. 

The Mach number evolution (Fig. [T) is similar in both codes 
and at all resolutions up to around 4 tj,, at which point all cal- 
culations show variations of order 5% from each other. No clear 
trends either between codes or with resolution are apparent, indi- 
cating that the variation observed is due to the stochastic nature of 
fully-developed turbulence producing different though statistically 
similar evolution (see Fig. [3] for the evidence for this in the den- 
sity field). It is these intermittent fluctuations that make turbulence 
code comparisons based on snapshot-to-snapshot comparison diffi- 
cult, because the instantaneous turbulence field will quickly diverge 
between different codes in the fully developed regime because of 
the chaotic nature of the turbulence (see projections and slices for 
t > 2td comparing the SPH and grid results). 

The evolution of maximum density (Fig.|2]( shows strong time 
variability in all six calculation s, similar to the resul ts shown in 
IKritsuk et al] feOOl Fig. 2) and lFederrath et al.l J2009I . Fig. 2). For 
isothermal flows arbitrarily large density fluctuations can be pro- 
duced, though with vanishingly small probability as can be in- 
ferred from the log-normal form of the PDF. This simply reflects 
the highly intermittent nature of the density fluctuations in super- 
sonic, turbulent flows. The maximum density is a clear function of 
resolution in each code, showing no signs of convergence, as one 
might expect seeing as we are sampling the very highest data point 
in the PDF. The results also demonstrate the evidently higher mass 
resolution in the SPH code: at 128 3 particles the maximum density 
resolvable in SPH is roughly similar to that resolved at 512 on 
the grid. Using 512 3 SPH particles the maximum density resolved 
at RMS Mach 10 is roughly three-and-a-half orders of magnitude 
above the mean density which one might therefore expect to be 
similar to the mass resolution in a 2048 3 grid-based calculation. 

3.2 Density field 

3.2.1 Projected density fields 

The projected column density fields at the highest resolution (512 3 ) 
are shown for PHANTOM, FLASH and the density field computed 
from the FLASH tracer particles in Fig. [3] (left, middle and right 
columns, respectively), at intervals of At = ltd for the first four 
dynamical times. The column density plots for SPH have been 
produced directly from the particles to a 2D pixel map using the 



SPLASH visualisation tool jPricdl2007h whilst the grid based re- 
sults have been integrated through the grid. We show the integration 
through the ^-direction in the codes. 

At early times the calculations show clear agreement in the 
location of individual shocks (t — ltd, top row) and in the devel- 
opment of large scale structures (t = 2td, second row). By t = 4td 
(fourth row) there is no longer clear correspondence even at the 
largest scales between codes, in agreement with the observed devi- 
ations in the time evolution of the RMS Mach number around this 
time (Fig.[TJ. 

In terms of resolution, high-density structures appear better 
resolved in the SPH calculations at the same number of computa- 
tional elements. However, the grid results tend to show better reso- 
lution of features in low density regions, as one might expect since 
in SPH the resolution is preferentially shifted away from low den- 
sity regions towards high density regions. 

The excellent agreement between codes in the development 
of individual shock structures within the first dynamical times (top 
row) enabled us to make a very detailed comparison of features be- 
tween codes which proved to be very helpful in the comparison pro- 
cess. In particular it highlighted that, with the parameters we were 
initially using, some of the dense structures created by the collision 
of one or more shocks were rapidly losing definition in the SPH 
results, resulting in a noisy density field that was rather unlike the 
grid results (this is shown in more detail in Appendix[A}- The prob- 
lem could be easily traced to be caused by particles penetrating or 
"overshooting" the shock front in these high Mach number shocks, 
a problem which the non-linear /3 (von-Neumann-Richtmyer) term 
in t he SPH artificial v iscosity (Eqs. 191 1 Ot was designed to prevent 
(see lMonaghan|[l989l) . The problem was thus easily fixed by us- 
ing a larger value for f3 V i SC . We have therefore used f3 V i SC = 4 
throughout the paper, rather than the nominal f3 V i SC — 2 which is 
widely used — and sufficient — for low Mach number calculations. 
It should be noted that this makes very little difference to the over- 
all dissipation rate since the linear viscosity term (a) dominates the 
numerical dissipation rate almost everywhere except at very strong 
divergence in the velocity field (where particle penetration can oc- 
cur). 

Comparing the projected column density fields calculated us- 
ing the tracer particles in the FLASH calculation (right column) to 
the grid-based density field (centre column), showing a wealth of 
sub-grid structures, suggests that the tracer particles have the abil- 
ity to provide a truly staggering improvement in resolution in the 
density field. The improved resolution is all the more remarkable 
considering that the tracer particles are merely advected with the 
grid-based velocity field at essentially no extra computational ex- 
pense. 

3.2.2 Cross section slices 

Column density plots such as those shown in Fig. [3] in general 
tend to highlight dense features, since all structures along the line 
of sight contribute to the projected field (which also tends to be 
the case in observations). The features in column density plots are 
therefore reflected by statistics such as the PDF (Figs. |6}{8} and 
quantities such as the maximum density (Fig. [2). However volu- 
metric quantities such as the volume filling factor of the material 
and the velocity field, reflected in statistics such as power spectra 
and structure functions, are better illustrated by cross-section slices. 
For this reason we show cross section slices of density at the mid- 
plane of the computational domain (z = 0.5), showing a resolution 
study of the initial shock development at 1 dynamical time (Fig. [4j( 
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Figure 3. Projected column density in the PHANTOM (SPH, left) and FLASH (grid, centre) calculations at a resolution of 512 3 particles/grid cells, showing the 
evolution over the first few dynamical times (top to bottom), together with the density computed from the tracer particle positions in the FLASH calculation 
using an SPH density estimate (right panel). After 1 dynamical time (top row), there is clear correspondence in individual shock structures between the SPH 
and the grid code, whilst after 2 dynamical times (second row) there are similar large scale features. However by 3 or 4 dynamical times (third and fourth 
row) only a weak correlation even between large scale features is observed. Dense features are in general better resolved in the SPH calculations at equivalent 
resolutions (in a number of particles=number of grid cells sense), whilst the grid-based calculations tend to better resolve features in low density regions (see 
also Fig.[5J. The increased resolution of sharp features in the tracer particle density fields (right column, compared to centre column) suggest a remarkable 
ability for the tracer particles to provide information on sub-grid scales at essentially zero additional computational cost. 
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Figure 4. Cross section slice of the density at the box midplane (z = 0.5) after 1 dynamical time, for three different resolutions (128 3 , 256 3 and 512 3 in grid 
cells/particles, top to bottom) using PHANTOM (SPH, left), FLASH (grid, centre) and for the density calculated from the tracer particles in the grid calculation 
using an SPH summation (right panels). 




Figure 5. Cross section slice of the density at the box midplane (z = 0.5), as in Fig.|4]but here shown after 10 dynamical times and showing only the highest 
resolution calculations (512 3 ). The FLASH calculations (grid, centre) shows better resolution in low density regions compared to PHANTOM (SPH, left). In 
evolved snapshots the tracer particles appear strongly clustered in high density regions and almost completely absent from the voids (right panel). 
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and a comparison of the evolved snapshots at the end of the simu- 
lations (t/td = 10), showing only the highest resolution (Fig. O. 
The plots show the density field using PHANTOM (left columns in 
Figs. [4] and [3), FLASH (centre column in both Figs.) and for the 
tracer particle density field computed from the FLASH calculations 
(right columns). 

Figs. E] and [5j show clearly that the grid results are better re- 
solved in low density regions. The resolution in the SPH calcula- 
tions is concentrated towards high density regions which fill rela- 
tively little of the volume. Comparing individual shock structures 
in Fig. [4] shows that in general the shocks have better definition 
in FLASH, with the shock widths in the highest resolution PHAN- 
TOM calculation similar to those obtained at 256 3 in FLASH. This 
is as might be expected given the relative crudeness of the shock 
capturing scheme (artificial viscosity) in the SPH code compared 
to the PPM shock capturing scheme JColella & Woodwar j[l984l) 
employed in FLASH. In the more evolved snapshots (Fig. |5_}, the 
grid results show many well-defined shock features in low density 
regions that are much less well resolved in the SPH calculations. 

Some numerical artefacts are visible in the lowest resolution 
SPH calculations in the earliest snapshot (t = ltd, top left panel of 
Fig. due to the "breaking" of the initial regular lattice on which 
the particles were placed as it is distorted by the flow. Interestingly 
similar artefacts are visible — and more accentuated — in the low 
resolution tracer particle plots (top right panel). These effects are 
not obviously visible either in the SPH or the tracer particles at 
higher resolution (middle and bottom rows of Fig. [4]l or at later 
times (Fig. [5} once the particles have adopted a more 'natural' ar- 
rangement. There are no obvious artefacts at low resolution in the 
grid based calculation. 

Density slices calculated from the tracer particles in the FLASH 
calculation are shown in the rightmost panels of Figs.|4]and|5] using 
the SPH density summation ((6} iterated self-consistently with the 
smoothing length according to 10. Comparison with the grid-based 
density field at 1 dynamical time, the tracer particles appear to sub- 
stantially increase the resolution in high density regions. Whilst 
most of the features have close correspondence to those visible in 
the grid slices (centre column of Fig. |4j, it is notable that a dense 
shock structure appears in the lower part of the t/td = 1 snap- 
shots at all resolutions that is completely absent from both the SPH 
and grid density fields. The absence of this feature even at 512 3 
in the centre column of Fig. [4] yet clearly present at 128 3 in the 
tracer particles, suggests that it may be an artefact of tracer parti- 
cles clustering below the grid scale. At later times (right panel of 
Fig. [5j this is even more evident by the fact that the tracer parti- 
cles are strongly concentrated in high density regions and largely 
evacuated from low density regions (i.e., large parts of the panels 
are saturated at the density floor of the plot due to the absence of 
a contribution from tracer particles even with iterated smoothing 
lengths). 

The difference between the density slices and column den- 
sity plots shows that in general, for similar numbers of compu- 
tational elements (not the same as equal computational expense), 
SPH codes are better at resolving dense structures (highlighted by 
projections through the volume), whilst grid codes are better at re- 
solving volumetric structures (highlighted by slices though the vol- 
ume). 



3.3 Probability Distribution Functions 

Many studies have demonstrated that the density probability 
distribution function (PDF) in supersonic turbulence is well 
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where the mean of the logarithm of density In p is related to the 
standard deviation a of In p by 



lnp = -a/2. 



(14) 



The appearance of a log-normal form in isothermal flows can 
be understood analytically as a consequence of the multi- 
plicative central limit theorem assuming that individual density 
pertu rb ations are independent and random JVazquez-Se madenil 



1994 IPassot & Vazquez-Semadenil 1 19981: iNordlund & PadoanI 
1999). In physical terms this means that density fluctua- 
tions at a given location are constructed by successive pas- 
sages of sh ocks with a jump amplitude in d ependent of the lo- 
cal density dBallesteros-Paredes etafl 120071 : iKritsuk et all 120071: 
IFederrath et alj|20ld) . Furthermore the width of the PDF for In p 
is found to be related to the rms Mach number according to 



ln(l + fc 2 M 2 l 



(15) 



where the factor b m 1/2 has been suggested by ea rly numerical 
experi ments (e.g. lPadoan et al.lll997h . More recently IKritsuk et al 

(2007) find a much lo wer value of b « 0.26 whilst Bee tz et al 

(2008) find b = 0.37. iBruntl d2010h has also recently measured 
b — 0.5 ± 0.05 in the Taurus Molecular Cloud, based on a method 
for inferring the 3D variance from 2D obs ervat ions developed 
by Brunt et al. J2010h . IFederrath et alj ( f2008l) and IFederrath et all 
( 1201 Of) reconcile these results by showing that the width of the 
PDF depends not only on the RMS Mach number but also on 
the relative degree of compressible and solenoidal modes in the 
forcing, with b — 0.33 appropriate for purel y solenoidal forc- 
ing an d b=l for purely compressive forcing. ILemaster & Stone] 
(2008) — performing calculations at a range of Mach numbers 
— suggest that the relationship l |15| > should be adjusted, finding 
a 2 = -0.72 In (l + 0.5M 2 ) + 0.20 from a three-parameter fit 
for hydrodynamic turbulence. For the purposes of the comparison 
at hand we simply fit the PDFs using a single parameter b based on 
Eqs. fUMB}. 



3.3. 1 Volume weighted PDFs 

Time-averaged PDFs of s = ln(p/po) for the three SPH calcula- 
tions and three grid calculations are shown in Figures [6] [7] and [8] 
The plots show the time average of individual PDFs computed at 
intervals of At = td/ 10, starting from 2td when turbulence is rea- 
sonably well established (see Figures [T] [3}. This gives a total of 81 
snapshots used in the averaging procedure. For reference we com- 
pute the PDF for each code with a bin width of 0. 1 in In p, with the 
first bin starting at (In p) m in = —12. 

For the grid results the volume weighted PDF is constructed 
simply by binning the grid cells according to the value of In p. To 
obtain a volume weighted PDF in SPH it is necessary to weight 
the contribution of each particle by the volume element associated 
with that particle, m / p, which for equal mass particles is simply in- 
versely proportional to the density. To construct the PDF we there- 
fore bin each particle according to the value of In p and add a con- 
tribution of 1/p to the bin, normalising the resultant PDF such that 
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Figure 6. Time-averaged Probability Distribution Function (PDF) of the 
logarithm of the density field s = In p from the PHANTOM (SPH, black 
lines other than dotted) and FLASH (grid, red lines) calculations, each at 
resolutions of 128 3 , 256 3 and 512 3 particles/grid cells. The PDFs are av- 
eraged over 81 snapshots evenly spaced between t/t^ = 2 and t/tj, = 10. 
Here we show the PDFs on a linear scale to highlight the change in posi- 
tion of the peak value as a function of resolution. We have also plotted, as 
dotted black lines, the best fit (to the peak) log-normal distributions from 
equations 1 1 3i and )15t using b = 0.5 (best fitting the 512 3 grid results) 
and b = 0.33 (best fitting the 512 3 SPH results). The SPH and grid results 
show complementary trends with resolution. 
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Figure 7. Probability Distribution Function (PDF) of In p, as in Fig. [6] but 
here shown on a logarithmic scale. The PDFs are log-normal to good ap- 
proximation for ~ 3 — 4 orders of magnitude in density either side of the 
mean, demonstrated by the best fit log-normal distributions (fitted to the 
peak in Fig. [6} given by the dotted black lines using b = 0.5 and b = 0.33 
to fit the 512 3 grid (solid red line) and 512 3 SPH (solid black line) results 
respectively. 



the integral over all bins (i.e, the total probability) is unity. This is 
different to the procedure used to construct density PDFs from SPH 
partic les used both in the Pot sdam comparison dKitsionas et all 
2009) (see howeve r Fig. 11 inlKitsionas et alj|2009 | ) and by pre- 



Klesse: 



by pre- 



vious authors (e.g. Vazquez-Semadeni et alj |2003l : [ 
iMac Low et al ] | 1998b . whereby the SPH results were first interpo- 
lated to a grid and a PDF constructed as above for the grid-based 
results. The main disadvantage to interpolating to the grid is that 
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Figure 8. Time averaged PDFs, as in Figures|6]and|2]but here showing the 
tails of the distributions at very low and very high densities. The SPH code 
(black lines other than dotted) resolves the PDF to much higher densities 
than the grid code (red lines), though the grid results correspondingly extend 
further into the low density regime. The log-normal distributions (dotted 
black lines) are no longer a good approximation to the distribution very far 
from the mean, in particular in the high density tail (where the 256 3 and 
512 3 SPH results appear to show convergence). 



the part of the high density tail of the SPH calculation that falls be- 
low the grid scale is removed. To retain this tail requires that PDFs 
should not be constructed from SPH particles by interpolating to a 
grid, though this is a perfectly valid procedure for computing volu- 
metric quantities such as power spectra (see £|3.41 l. 

The PDFs thus constructed (Fig.O show clearly a log-normal 
distribution in agreement with many previous calculations and with 
theoretical expectations (see above). Whilst the results are broadly 
similar for all calculations in the central regions around the mean 
density and in overall shape, clear differences may also be observed 
between codes and with resolution, particularly in the tails of the 
distribution (Fig.[8](. At the low density end (Fig.(8]( the grid results 
tend to show a wide, low density tail, with probability densities 
decreasing with resolution. By contrast, the SPH results show a 
narrower low density tail with probability densities that increase 
with resolution. 

Whilst both codes appear to be converging towards each other, 
it is clear that neither is well converged at the low density end 
(p/po < 0.01) at least for the resolutions used in this paper. Thus 
the low density tail should not be used to fit the PDF width from 
either grid or SPH codes alone at these resolutions. Instead we mea- 
sure the PDF width using the best fit around the mean value (i.e., 
based on Fig. [6]». The resulting best fit log-normal distributions are 
plotted as dashed lines in Fig. [6] corresponding to b = 0.33 for the 
512 3 SPH results and b = 0.5 for the 512 3 grid results. As it turns 
out, the best fitting log-normal distributions also provide a good fit 
to the low density tails (Figures [7] and [8}, though the same is not 
true at very high density (see Fig.[8jl. 

At the high density end the trend with resolution for both 
codes is for the PDF in the high density tail to increase slightly, 
increasing the PDF width. The FLASH results show a stronger trend 
with resolution at high densities and appear to converge towards the 
PHANTOM results, with rough equivalence between the 512 3 grid 
based results and the 128 3 SPH results in resolving the high den- 
sity tail, similar to what is observed for the evolution of maximum 
density in Fig. [2] Interestingly, the SPH results appear converged 
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between the 256 3 and 512 3 calculations for densities up to around 
p/po ~ 10 3 (though not to the best fit log-normal, see below), with 
the primary effect of the additional resolution at 512 3 being to ex- 
tend the tail to lower probabilities (but with similar overall width). 

That the SPH results appear close to converged between the 
256 3 and 512 3 and that the grid results are converging towards 
them suggests that a value of b ~ 0.35 — 0.4 would be the con- 
verged value. This is similar to the results of iBeetz et al. I d2008h 
for solenoidal forcing (6 = 0.37), and close to the prediction of 
6=1/3 for solenoidal forcing from the heur istic model for the b 
param eter presented in Federrat h et al.l d2008l) and lFederrath et ail 

hold) . 

The reason for the discrepancy at high density merits some 
consideration given the importance of this regime in relation to star 
formation. The most straightforward conclusion is that the conver- 
gence we find is merely incidental and that performing calculations 
at higher resolution would resolve the remaining discrepancy. Cer- 
tainly, resolution requirements on the PDF become greater at higher 
Mach numbers because of the stronger time-variability in the tails 
of the distribution. This is evident from the strong fluctuations in 
pmax (Fig. [2} that are up to an orde r of magnitude larg er than 
the fluctuations shown in Figure 2 of iKritsuk et alj d2007l) . There 
may also b e differences due to the random forcing algorithm. In 
particular Federrath et al. (2010), using the same (solenoidal) forc- 
ing algorithm employed here, also find a small deviation from log- 
normality at the high density end, though smaller than we find since 
they employ 1024 3 elements and a lower Mach number (Mach 
6). They also found strong deviations from log-normality when a 
compressible forcing was applied, indicating that the PDF at high 
densit ies is quite sensitive to the forcing employed. Federra th et al.l 
(2010) discuss intermittency as a cause of non-Gaussian PDF 
tails. They quantify and discuss the Mach num ber-density corre- 
lations as the key to non-Gaussian PDFs (see Vazquez-Semadeni 
Il994l: IPassot & Vazquez-Semadenill 19981) . iFalgarone et alj Jl994h 
and iHilv-Blant et alj {2008) also find strong intermittent fluctua- 
tions in their molecular cloud observations and attribute them to 
a fundamental proper ty of turbulence, i.e., intermittency (see, e.g., 
IShe & Levequell 1994b . 

In principle it is also possible to compute the PDFs from the 
tracer particle density field calculated with the SPH summation. 
However we find that the resulting PDFs show a strong deviation 
from a log-normal distribution, particularly in the high density tail 
(much stronger than those seen in Fig. [8] and in the opposite direc- 
tion), due to the manner in which tracer particles tend to cluster 
in high density regions at later times (see discussion in i)3.2.2| and 
right panels in Fig. [5}. 



3.4 Power spectra 

Padoan & Nordlunc] J2002h derive a relationship between the mass 
distribution of dense cores to the slope of the kinetic energy power 
spectrum in supersonic, super-Alfvenic turbulence via the relation 



N (m)d log m <x m ^dlogm, 



(16) 



where /3 is the slope of the kinetic energy power spectrum assumed 
to be a power law of the form 

E{k) oc k' 13 , (17) 

where /3 would be ~ 5/3 according to the iKolmogorovl J 1 94 lh 
phenomenology for incompressible turbulence. Supersonic turbu- 
lence is generally found to have a power spectrum closer to the 
completely pressure-free shock-dominated turbulence produced by 




Figure 9. Velocity power spectra (^V 



), shown compensated by k 2 as an 



average over 81 snapshots evenly spaced between t/t^ = 2 and t/t^ = 10 
for the SPH (solid, black) and grid (dashed, red) calculations at the three 
different resolutions, as indicated. For calculations at or above 256 3 compu- 
tational elements in either code the results in the scaling range 7 < k < 11 
are consistent with a slope slightly shallower than Burgers' value of k~ 2 , 
between A; -1 93 and k~ 198 depending on whether or not the 7 < k < 11 
region is interpreted as inertia] range or a bottleneck effect. 



solving Burgers' equation, f3 — 2, implying dominance of the non- 
linear advection term (v ■ X7)v over the pressure gradie nt in the 
equati on of motion. A value consistent with the observed Salpeter 
d 19551) slope for the IM F of N(m)d\og m oc m~ " 4 / 3 dlo K m re- 
quires P w 1.75, which IPadoan & Nordlunj d2002l) suggest arises 
from supersonic MHD turbul e nce in the super-Alfvenic regime. 

Recently IKritsuk et alj d2007l) have suggested that Kol- 
mogorov scaling may be applicable for highly compressible tur- 
bulence by assuming that the mean volume energy transfer rate 
pv 2 v/l, is constant, implying that E(k) oc k~ 5 ^ 3 holds for the 
variable w = (p 1 ^ 3 v) rather than simply for the velocity. We thus 
consider power spectra of both quantities. 

We have computed the power spectra for each code directly 
from gridded data using the same analysis script for both codes. 
For the SPH code, this means that the results have first been interpo- 
lated to grids of size 256 3 , 512 3 and 512 3 cells for the 128 3 , 256 3 , 
512 3 particle calculations respectively. The interpolation has been 
performed using a routine from t he SPLASH visualisation code, the 
details of which are described in |Priceld2007l) . The main disadvan- 
tage of interpolating SPH data to a grid is that any resolution in 
the SPH code on scales smaller than the grid scale is lost. This 
is most obvious when comparing quantities such as the maximum 
density on the particle data compared to the maximum density on 
the interpolated grid. For example the maximum density interpo- 
lated onto a 512 3 grid (max(p/p ) ~ 3 — 4 x 10 2 ) is a factor 
of ~ 3 — 4 lower than the maximum density from the particles 
(max(p/p ) ~ 1 - 2 X 10 3 ) for snapshots from the 512 3 SPH 
calculation, which would remove some of the information from the 
high density tail of the PDF as discussed above. Whether or not 
the interpolation procedure affects the power spectrum calculation 
can be determined simply by comparing the results from interpola- 
tions to different sized grids. We find that for the power spectrum, 
as one might expect given that it is a volumetric measurement, the 
power spectra are identical for different grid sizes apart from fc's 
very close to the grid scale. 
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Figure 10. Density- weighted velocity power spectra, as in Fig. [9] but for 
the quantity (p 1 ' 3 v) instead of the velocity field. There is tentative ev- 
idence of a flat portion in the compensated spectra at resolutions above 
256 3 in either code, suggesting a small scaling range in the region of 
4 S; k < 6 (dotted line) consiste nt with a Kolmogorov -like fc -5 / 3 scaling 
for this quantity as suggested by Kritsuk et al. 1 2007). This also suggests 
that 512 3 particles/grid cells is the minimum resolution requirement to de- 
termin e the power spectru m slope in the inertial range in either code (see 
also lFederrath et alj201ol Fig. C.l). 




Figure 11. Power spectra, as in Figs. l9land l 101 but here for the density field, 
computed as the spectrum of the density fluctuations 8p = p — p. Clearly 
the SPH code provides increased resolution in the density field at high k 
compared to the grid-based code for the same number of computational el- 
ements. The spectra have been compensated by k 2 , though it is evident that 
the slope of the density spectrum is not well described by a single power- 
law. 



Power spectra 

£(k) = iw(k) 2 , (18) 

for an arbitrary vector field w are constructed from the 3D Fourier 
transform, 



w(k) 



(2tt) 3 



w(x)e d x, 



(19) 



dard procedure of summing E( k) in bins according to the norm of 
the wavenumber, |fc| (as in e.g. lKitsionas et alj2009l) . 



3.4.1 Volume weighted velocity power spectra 

The velocity power spectrum (that is, where w = v), averaged over 
81 snapshots as described above for the PDFs, is shown — com- 
pensated by k' 2 — in Fig. [9] for the three SPH calculations (solid 
black lines) and for the three grid calculations (dashed red lines). 
The results, similar to previous studies, show a peak at the driving 
scale k ~ 2, a power law slope (flat in these compensated spectra) 
between k ~ 4 and up to k ~ 12 in the highest resolution calcu- 
lations and an extended dissipative tail at large k. Both codes show 
a power law slope close to the pressure-free Burgers model slope 
of k~ 2 . The trend with resolution is for both codes to tend towards 
a slightly shallower slope, with the highest resolution calculations 
consistent with ~ fc -1 ' 95 — fc -1 ' 98 for both codes at 512 3 , depend- 
ing on whether or not the results for 7 < k < 12 are interpreted as 
inertial range or as a bottleneck effect. The extent of the power-law 
scaling range of the velocity power spectrum appears to be very 
similar in both codes at equivalent (number of grid cells = num- 
ber of particles) resolutions, though the SPH results show a faster 
drop-off tow ards higher k in the di ssipative tail in agreement with 
the results bv lKitsionas et al.l d2009h . 

These results are en tirely consistent with the power law slopes 
obtained in Figure 2 of Pado an et al. I d2007l) . However they dif- 
fer strongly from the low resolution SPH an d TVD results shown 
in Fig. 2 of iBallesteros-Paredes et al" I d2006l) . shown compensated 
in Fig. 8 of iPadoan et al.l ( 120071) . In particular our SPH results 
show power law slopes consistent with j3 — 2 or shallower, 
which is in stark contrast to t he f3 = 2.7 — 2.9 obtained by 



The angle-averaged power spectrum was then obtained by the stan- 



iBallesteros-Paredes et alj ((2006), albeit at lower Mach numbers (6 
and 3 respectively). This can be understood primarily as an effect 
of the numerical resolution, since we observe a significantly steeper 
slope in our low resolution SPH calculations shown in Fig. [9] As 
noted in the introduction however, the resolution in our 'low' reso- 
lution SPH calculation, employing 2.1 million particles, is already 
an order of magnitude higher than tha n the ~ 200, 000 parti- 
cles u sed for most of the calculations inlBallesteros-Pared es et all 
(2006) (they perform one 'high' resolution calculation using « 10 
million particles, the power spectrum of which is not shown in their 
paper, though it is used to derive a core mass distribution). We also 
find excellent agreement b etween both the SPH and grid-ba sed re- 
sults which is in contrast to Ballesteros-Par edes et al.l 12006) where 
the measured slopes differ between their codes which had been 
attributed to the very different properties of the SPH and TVD 
schemes used in their paper. 



3.4.2 Density weighted velocity power spectra ( p 1 ^ 3 v) 

Fig. QO] shows the time-averaged power spectrum, computed as 
above, of t he quantity w = p 1//3 v, which according to the hy- 
pothesis of iKritsuk et alj J2007h is the quantity that for supersonic 
turbulence should show a Kolmogorov-like scaling of fc~ 5 / 3 . We 
have therefore compensated the spectra by fc 5//3 in order to assess 
whether or not this can be supported on the basis of our calcula- 
tions. The spectra in Fig.[l0]show similar generic features to those 
observed in Fig. [9] 

The lower resolution calculations (128 3 and 256 3 ) in both 
codes show slopes that appear shallower than f3 = 5/3, though 
the convergence of both codes is towards a steeper slope with res- 
olution. In particular the 512 3 calculations with both codes show a 
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small, flat region in the compensated spectrum between k ~ 4 and 
k ~ 6 that may be interpreted as a resolved inertial range ( shown 
by the dotted black line) consistent with the iKritsuk et aP d2007h 
scaling, with a 'bottleneck effect' for k > 6 extending into the 
dissipa tive tail. For the grid cod e, this is consist e nt wi th the find- 
ings of IKritsuk et"ai] J2007t) and lFederrath et al.1 d201dh that 512 3 
is roughly the minimum resolution required to resolve the inertial 
range, and we conclude that a similar requirement holds for SPH. 

Interestingly, the difference between codes in rate of drop off 
in the high-fc dissipative tail for the p 1//3 v spectrum is the reverse 
of what occurs for the velocity field alone (comparing to Fig. [9^. 
That is, for p 1 ' 3 v the grid-based results drop off much faster at 
high-fc than the SPH code, whereas the SPH code drops off faster 
in the velocity spectrum. We interpret this as being due to the fact 
that SPH has better mass resolution for equal numbers of compu- 
tational elements (reflected in the density field) but worse volume 
resolution (reflected in the velocity field). This is further evident in 
the spectrum of the density field alone, discussed below. 



in the computational domain, with a constraint to achieve similar 
sampling of all lags. For the structure function calculations we typ- 
ically use ~ 250, 000 sample points. For the SPH results we have 
firstly computed structure functions by interpolating the SPH data 
to a grid, as described above for the power spectra. Whilst this en- 
sures that the analysis step is identical to the grid-based results, this 
procedure has the disadvantage that the SPH densities and veloci- 
ties are smoothed by the interpolation step. For comparison we have 
therefore computed structure functions also directly from the par- 
ticle data. In the particle version we iterate the structure functions 
to completion by progressively increasing the number of sample 
points by a factor of 2 and check that the error between the final 
result and the structure functions using half the number of sample 
points is small (typically we require the RMS error on the highest 
order structure function to be < 1%). In practise this produces a 
similar sample size to that used for the grids i.e., a few hundred 
thousand sample points for a 512 3 calculation, but can be more ef- 
ficient if the velocity field is well ordered (e.g. at early times). 



3.4.3 Density power spectra 

Fig. QTJ shows the power spectra of the density field, computed as 
the spectrum of density fluctuations Sp = p — p such that the in- 
tegral under the power spectrum gives the density variance. Whilst 
not well fit by a single power law we show the spectra compensated 
by k 2 to facilitate the comparison. The overall shape of the spectra 
is similar between codes on large scales (small k), though it is clear 
that the SPH code shows a much higher resolution in the density 
field, falling much more slowly at high k compared to the grid re- 
sults at comparable resolutions and explaining the results in Fig. 1101 
as intermediate between Figures [9] and QTJ The SPH density spec- 
trum at 128 3 particles appears better resolved than the 256 3 grid 
spectrum though not as well resolved as the 512 3 grid spectrum, 
lying somewhere between the two. 

3.5 Structure Functions 

We have measured the structure functions, 

S p (0 = (!w(r + /)-w(r)| p ) (20) 

from our calculations, where in this paper we show orders up to 5, 
i.e., p — 1, 2, 3, 4 and 5 and the variable w is either the velocity 
w = v or the quantity w = (p 1//3 v) as for the computation of 
the power spectra in the previous section. We compute intrinsically 
three dimensional structure functions (rather than the approximate 
'directionally split' versions used e.g. in the KITP comparison) by 
sampling pairs of either grid cells or SPH particles, computing the 
relative velocity difference and adding the result (the absolute value 
to the power p) to the appropriate bin corresponding to the spatial 
separation of the points (the "lag" — given here in computational 
units i.e. in terms of the box length L). Dividing by the total number 
of contributions to each bin produces the average, as indicated by 
the angle brackets (...). We have computed the structure functions 
for both the transverse (w perpendicular to the line 1 joining the 
pair of points) and longitudinal (w parallel to 1) components of w 
in each case. 

Calculation of the structure functions is computationally ex- 
pensive, in principle requiring a summation over 0(N 2 ) pairs of 
points. A more efficient calculation can be made by selecting only 
a subsample of points. We achieve this by randomly selecting a 
fixed number of points (either grid cells or particles) and comput- 
ing the pair-wise interaction of these sample points with all points 



3.5.1 Velocity structure functions 

The structure functions for the velocity field (w = v) at all com- 
puted orders (p = 1..5) are shown in the left hand side of Fig. 1121 
for the 512 3 calculations with both codes. The structure functions 
show the same general features in both codes, showing slopes that 
appear power-law like in the range I ~ 0.03 to I ~ 0.3 (though see 
below), a flat structure around the driving range (0.3 < I < 0.5) 
and a fall-off in the dissipative tail (/ < 0.03). The transverse struc- 
ture functions appear flatter around the driving scale compared to 
the longitudinal version, consistent with the fact that we have used 
purely solenoidal forcing. Differences between codes mainly ap- 
pear in the dissipative tail, as observed also in the power spectra 
(Fi gures |9l and 1 1 01) . In the velocity structure functions, as with the 
velocity power spectrum, it appears that the SPH results fall off 
faster in the dissipative tail (comparing solid and dot-dashed lines 
in the left hand part of Fig. 112b . though the results computed di- 
rectly from the particles (dashed lines) appear to show a slower 
fall-off, indicating that at least some of this may be an artefact of 
the interpolation procedure. However the structure functions com- 
puted directly from the particles also show much shallower slopes 
at the lowest orders (p = 1..4), so it is not clear that the structure 
functions computed in this manner are truly comparable. 



3.5.2 Density weighted velocity structure functions ( p 1//3 v) 

Following IKritsuk et al] d2007l) , we have also computed structure 
functions for w = (p 1/,3 v). These are shown on the right hand side 
of Fig. [12] showing — given the above — only the versions com- 
puted identically from a grid (i.e., interpolated in the SPH case). 
Again, the qualitative features are in good agreement between the 
two codes — in particular the structure functions for p 1//3 v show 
a power-law scaling range over a much wider range of distances, 
0.01 < I < 0.3. As was the case for the power spectra (see Fig.UOIl, 
with the density-weighting the grid code (solid lines) shows a faster 
fall-off in the dissipative tail compared to the SPH code (dashed 
lines), in contrast to the pure velocity version (comparing the right 
hand and left hand panels in Fig.ll2t. 
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Figure 12. Structure functions S p = (|w(r + /) — w(r)| p ) as a function of spatial separation (the lag, I), for orders p = 1,2,3,4 and 5, showing the time 
average over 81 snapshots evenly spaced from = 2 to = 10 in each case. Structure functions are shown for the longitudinal (w || 1, top) and transverse 
(w _L 1, bottom) components of the velocity field, w = v (left) and for the mass-weighted velocity w = p 1 / 3 v (right). The results from the FLASH (grid) 
512 3 calculation are given by the dashed red lines whilst the PHANTOM (SPH) results at 512 3 are shown by the solid black lines, where these have first been 
interpolated to a 512 3 grid in order to be analysed identically to the grid results. Results for the velocity field computed directly from the SPH particles for the 
512 3 calculation are plotted on the left figure using dashed black lines. 




Figure 13. Time-average of the measured structure function slope £ p , where S p oc I^p for the transverse velocity (left) and p 1 ' 3 v (right), averaged from 81 
snapshots evenly spaced between t/tj = 2 and t/i<j = 10, and plotted as a function of the structure function order p. The solid black lines correspond to 
the PHANTOM (SPH) calculations whilst the red dashed lines correspond to the FLASH (grid) calculations, each shown at 3 different resolutions (the slopes 
decrease with resolution in both codes). Error bars show the 1<t errors in the time-average of the measured slope. The corresponding time average of the 
structure functions themselves are shown in the lower panels of Fig. [T2l 
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3.6 Structure function scaling 

A more quantitative analysis can be made by performing a least 
squares fit to determine the structure function slope Q p assuming 
S p oc l^ p for each order p. The results of such a fit to the time- 
averaged structure functions shown in Fig.[T2] — performed over a 
fixed Al in the range 0.1 < I < 0.15 (see below) — are shown in 
Fig. [73] Plotted are the fitted slope ( p at each order as a function 
of p, for both the velocity structure functions (left panel) and the 
density- weighted versions (right panel). The results from the SPH 
code, computed from the grid-interpolated structure functions, are 
shown in black, whilst the FLASH (grid) results are shown in red. In 
order to quantify the time variability in the value of the slope, we 
have also computed the standard deviation in the structure function 
slope at each order computed for each of the 8 1 snapshots discussed 
above. The 1 — a deviations are plotted as the error bars in Fig. 1 1 3 1 

For the velocity structure functions (left panel of Fig.|13t, rea- 
sonable agreement between codes is seemingly obtained for the two 
lowest order structure functions (p = 1 and 2) at the highest res- 
olution employed (512 3 , though see below). The grid code (red, 
dashed lines) appears close to (though not fully) self-converged at 
these orders (but not for p > 2) whilst the SPH code (solid, black 
lines) appears only close to convergence at the lowest order. The sit- 
uation appears slightly better for the p 1/ ^ 3 -weighted structure func- 
tions (right panel of Fig.[T3l. where the grid code (red, dashed lines) 
shows apparent convergence in the slopes up to p — 3 between the 
256 3 and 512 3 calculations. The SPH results (black, solid lines) 
appear to show only a small change between the measured slopes 
at all orders between the 256 3 and 512 3 results, though given the 
remaining disagreement between codes even at 512 3 it is not clear 
that true convergence has been reached in either. 

The somewhat large caveat to the results shown in Fig. 1 131 is 
our finding that the determined slope depends strongly on the range 
of scales (Al = l max — Imin) used to perform the least squares fit. 
This is discussed in more detail in Appendix iBl and demonstrated 
by Fig. lBll that shows the dependence of £ p on the employed Al for 
the transverse structure functions. It is clear from the fact that the 
measured slope steepens as Al is increased that relatively little true 
power-law scaling range is apparent in either code at the resolutions 
employed. This is particularly true for the velocity structure func- 
tions (left panel of Fig. IB It where either a very short (Al < 0.02) 
scaling range exists or none at all. The situation is marginally im- 
proved for the density-weighted structure functions (right panel of 
Fig. IB It . where the measured slope is less dependent on the fit- 
ting range (also apparent from Fig.|13t. However there remains a 
significant dependence on Al even for the lowest order p — 1. Fur- 
thermore, it is clear that any nominal convergence between 256 3 
and 512 3 results in individual codes seen in Fig.[T3]o;?Zy occurs if a 
relatively short range in / (Al < 0.05, i.e. 0.1 < I < 0.15) is used 
to perform the fit (and, as noted above, the measured slopes in this 
range also do not agree between codes). It is therefore difficult to 
assert that either code produces a resolved scaling in either the ve- 
locity or density weighted structure functions at any order at these 
resolutions. 

Given the above, we have not attempted a detailed compari- 
son with phenomenological m odels for the structure f unction scal- 
ing, such as those proposed by She & Leveque ( 1994) (revised for 
supersonic turbulence by lBoldvrevll2002l and lSchmidt e t al. 2008). 
Nevertheless, although the true structure function scaling remains 
to be determined, it is cl ear that the scalin g for both v and 
does not match either the Kolmogorov d 1 94 lh (£ p oc p/3) or Burg- 
ers (fp = 1 for p > 1) expectations, and also does not match the 



Boldyrevl J2002l) theory, though the relative scaling we have mea- 
su red does compare re asonably well to the parameters suggested 
bv lSchmidt etall J2008I) . 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have performed a detailed comparison of the statis- 
tics of supersonic turbulence using high resolution numerical simu- 
lations with two fundamentally different codes: FLASH, an Eulerian 
grid-based code and PHANTOM, a Lagrangian particle-based SPH 
code. Despite the very different numerical methods we find in gen- 
eral very good agreement between the codes on the many aspects 
of supersonic turbulence, though it is clear that neither code shows 
results that are fully converged at 512 except for a small scaling 
range in the velocity power spectrum with both codes and some in- 
dication of convergence near the peak and at the high density end 
of the PDF in the SPH case. 

We find good agreement in the fact that hydrodynamic turbu- 
lence at Mach 10 has a velocity power spectrum with a slope around 
E(k) oc fc -1,95 , very close to Burger's value of E(k) oc k~ 2 , 
confirming many previous results using only grid-based codes (e.g. 
iPadoan et al.l2007l ; lFederrath et alj201fj|) . At the highest resolution 
employed, both codes support the idea that the power spectrum of 
the mixed quantity p 1,/3 v shows a Kol mogorov-like scaling in the 
power spectrum of k 5 ^ 3 , as proposed bv lKritsuk et al J 12007). Both 
codes agree that the PDF of the logarithm of the density shows a 
distribution that is log-normal to good approximation for densities 
around 3 — 4 orders of magnitude either side of the mean density, 
with the factor b relating the width to the Mach number measured 
to lie between b — 0.33 to b — 0.5, with the converged value ex- 
pected to lie somewhere around b « 0.4. However, both codes also 
show that there are significant deviatio ns from a log-normal distri- 
bution at high density, as also found bv lFederrath et al.l d2010h . 

Our conclusions regarding resolution and computational re- 
quirements are as follows: 

(i) For measuring volumetric statistics such as the power spec- 
trum slope and structure function scaling, we find that SPH and 
grid codes give roughly comparable results when the number of 
SPH particles is approximately equal to the number of grid cells. 
This means that SPH codes are not well suited to measuring such 
quantities since for this kind of problem an SPH code will be sig- 
nificantly more expensive than a uniform-grid implementation at 
similar numbers of computational elements (we find the cost of the 
SPH calculations with PHANTOM similar to the cost of running the 
grid code (FLASH) at twice the resolution, i.e., with 256 3 particles 
« 512 3 grid cells in terms of CPU time). 

(ii) On the other hand the SPH code was found to be better at 
resolving dense structures, giving maximum densities at a resolu- 
tion of 128 3 particles that were similar to the maximum densities 
resolved in the grid code at 512 3 cells, reflected also in the high 
density tail of the PDF. SPH is therefore a more efficient method in 
this regard, which is an important reason why it is frequently used 
for studying star formation (for which a grid based code requires 
adaptive mesh refinement, adding a significant computational over- 
head). 

(iii) At comparable resolutions (N par t ~ N ce u s ), SPH 
(PHANTOM) appears to be more dissipative (steeper dissipative tails 
in the power spectrum and structure functions) in the velocity field, 
but the reverse is true for the statistics of the density-weighted 
quantity p 1,/3 v where the grid (FLASH) results appear more dis- 
sipative. We attribute the former to the greater sophistication of the 
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shock capturing scheme in the grid-based code, and the latter to 
the better resolved power in the density field provided by the SPH 
code. 

(iv) The absolute values of the structure function slopes are not 
converged in either grid or SPH at 512 3 , requiring much higher 
resolution in order to make a meaningful comparison with scaling 
models. 

(v) In order to accurately simulate supersonic turbulence in SPH 
it is important to ensure that sufficient /3-viscosity is applied to pre- 
vent particle interpenetration in shocks. At Mach 10 we require 
Pvisc = 4 instead of the usual f3 V i SC — 2. 

(vi) We find that calculation of the sub-grid density field from 
the tracer particle distribution using an SPH summation can signif- 
icantly enhance the resolution of high density structures from the 
grid-based results. However it is unclear whether or not the statis- 
tics of these sub-grid structures can be used in a meaningful way, 
because of the manner in which tracer particles tend to cluster in 
high density regions. 

The above conclusions mean that the differences b e tween the 
SPH and grid based results discussed by Pado an et alj j2007h . at 
least in terms of the power spectrum slope, can be understood as 
a consequence of the low resolution employed in the SPH cal- 
culations rather than being due to an intrinsic difficulty with the 
method. Given this, it is also lik ely that the conclusions drawn by 
iBallesteros-Paredes et alj {2006) regarding the presence or other- 
wise of an emergent power law slope in the core mass distribu- 
tion should also be treated with caution. As a follow-up to this 
paper it would be interesting to examine the properties of dense 
clumps (or 'cores') in our cal culations using a clump-finding algo- 
rifhm in a similar m anner to IBallesteros-Paredes et alj d2006t) and 
|Padoanetal.l ( l2007h to see whether or not these differences can also 
be reconciled by improved resolution in the SPH results. 
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APPENDIX A: EFFECT OF /3 VISCOSITY IN THE SPH 
CALCULATIONS 

Fig. IA1 1 shows the results of three 128 3 particle SPH calculations 
after 2 dynamical times using /3 V i SC — 1 (top), /3 v i sc = 2 (mid- 
dle) and /3 V i SC = 4 (bottom) in the SPH artificial viscosity term 
d9t-d 10>. At this high Mach number using /3 V i SC < 4 means that 
particle interpenetration can occur at the shock fronts (i.e., particles 
overshoot the shock), resulting in excess noise in the density field 
compared to the grid-based results and compared to the SPH results 
at higher /3 v i 3C . In particular the shocks appear much sharper and 
well-defined at higher /3 V i SC , somewhat counter intuitively since we 
are adding more viscosity. Using /3 V i sc = 4 shows good agreement 
with the grid-based shock structures (c.f. Figs.[3]and[4]l. 



APPENDIX B: EFFECT OF FITTING RANGE ON THE 
MEASURED STRUCTURE FUNCTION SLOPE 

The best-fit structure function slopes, £ p , computed from a least- 
squares fit to the transverse velocity and density-weighted structure 
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Figure Al. Effect of beta viscosity in the SPH calculations. From left to 
right: column density at t/tj, = 2 in 128 3 SPH calculations using fi v i sc = 
1, fivisc = 2 and f) v i ac = 4. With /3 v i sc < 2 particle penetration occurs 
in the shocks at these high Mach numbers, causing them to lose definition. 
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Figure Bl. The dependence of the measured structure function slope £ p on 
the range in scales (lag) over which the least squares fit to the time aver- 
aged structure function is performed, where Al = lmax — Imin ana we 
have used a fixed l max = 0.15. The measured slope, shown for the three 
lowest-order structure functions (p = 1,2, 3) shows a strong dependence 
on the adopted fitting range, particularly for the velocity structure functions 
(left panel) and to a lesser extent for the density-weighted (p 1 / 3 v) versions 
(right panel). Thus, although the slopes appear to converge in individual 
codes for a fixed Al, there remains a strong dependence on Al, indicating 
that at a scaling range of constant f p is not resolved in either code at the 
resolutions employed. 



functions shown in Fig. [13] are shown in Fig. IB II as a function of 
the range in the lag Al = lmax — Imin used to perform the fitting. 
To produce the Figure we initially adopted a fixed lmax = 0.15, 
based on visual inspection of Fig. [T2] for the transverse case, and 
performed the fit down to l m in = lmax — A/, plotting the result- 
ing slope as a function of Al. The results show a strong steepen- 
ing of the measured slope as the fitting range is increased, indi- 
cating that the slope is not well-fitted by a single power-law. The 
p 1//3 -weighted structure functions (right panel) show a weaker de- 
pendence on Al than for the pure velocity equivalents (left panel), 
though the results from the two codes even in this case do not show 
convergence with each other. Convergence in each individual code 
is only obtained if a relatively short fitting range Al < 0.05 is 
adopted, but it is clear that the converged value will nevertheless 
change if a different Al is used. Neither can it be asserted that the 
adopted Al should be changed with resolution, since it is not clear 
that any scaling range has been resolved even at the highest reso- 
lution adopted in either code (making the "correct" choice Al — 
for the numerical resolutions tested in this code comparison). 



